
program atmos_model

!-----------------------------------------------------------------------
!
!  Main program for running a stand-alone atmospheric dynamical core.
!
!-----------------------------------------------------------------------

#ifdef INTERNAL_FILE_NML
use mpp_mod, only: input_nml_file
#else
use fms_mod, only: open_namelist_file
#endif

use   atmosphere_mod, only: atmosphere_init, atmosphere_end, atmosphere, atmosphere_domain

use time_manager_mod, only: time_type, set_time, get_time,  &
                            operator(+), operator (<), operator (>), &
                            operator (/=), operator (/), operator (*), &
                            set_calendar_type, set_date, get_date, days_in_year, days_in_month, &
                            NO_CALENDAR, THIRTY_DAY_MONTHS, NOLEAP, JULIAN, GREGORIAN, INVALID_CALENDAR

use          fms_mod, only: file_exist, check_nml_error,                &
                            error_mesg, FATAL, WARNING,                 &
                            mpp_pe, mpp_root_pe, fms_init, fms_end,     &
                            stdlog, stdout, write_version_number,       &
                            open_restart_file, lowercase,               &
                            mpp_clock_id, mpp_clock_begin,              &
                            mpp_clock_end, CLOCK_COMPONENT, set_domain, nullify_domain
use       fms_io_mod, only: fms_io_exit

use  mpp_mod,         only: mpp_set_current_pelist
use  mpp_domains_mod, only: domain2d
use       mpp_io_mod, only: mpp_open, mpp_close, MPP_ASCII, MPP_OVERWR, &
                            MPP_SEQUENTIAL, MPP_SINGLE, MPP_RDONLY, MPP_DELETE

use diag_manager_mod, only: diag_manager_init, diag_manager_end, get_base_date

use  field_manager_mod, only: MODEL_ATMOS
use tracer_manager_mod, only: register_tracers
use       memutils_mod, only: print_memuse_stats
use   constants_mod,    only: SECONDS_PER_HOUR,  SECONDS_PER_MINUTE

implicit none

!-----------------------------------------------------------------------

character(len=128), parameter :: version = &
'$Id$'

character(len=128), parameter :: tag = &
'$Name$'

!-----------------------------------------------------------------------
!       ----- model time -----

   integer :: calendartype = INVALID_CALENDAR
   type (time_type) :: Time, Time_init, Time_end, Time_step_atmos
   integer :: num_atmos_calls, na
   type (time_type) :: Time_tmp ! used to facilitate some operations on time_type variables.
   integer :: days_tmp          ! used together with Time_tmp 

! ----- model initial date -----

   integer :: date_init(6)

! ----- timing flags -----

   integer :: id_init, id_loop, id_end
   integer, parameter :: timing_level = 1

!-----------------------------------------------------------------------
   character(len=80) :: text
!-----------------------------------------------------------------------
   type(domain2d), save :: atmos_domain  ! This variable must be treated as read-only
!-----------------------------------------------------------------------

      character(len=32) :: calendar="no_calendar" ! valid options are "no_calendar", "thirty_day_months", "noleap", "julian", "gregorian"
      integer, dimension(4) :: current_time = (/ 0, 0, 0, 0 /)
      integer :: years=0, months=0, days=0, hours=0, minutes=0, seconds=0
      integer :: dt_atmos = 0
      integer :: memuse_interval = 72
      integer :: atmos_nthreads = 1

      namelist /main_nml/ calendar, current_time, dt_atmos,  &
                          years, months, days, hours, minutes, seconds, memuse_interval, atmos_nthreads

!#######################################################################

 call fms_init ( )
 call atmos_model_init 

!   ------ atmosphere integration loop -------

    call mpp_clock_begin (id_loop)

    do na = 1, num_atmos_calls

       call atmosphere (Time)

       Time = Time + Time_step_atmos

       if(modulo(na,memuse_interval) == 0) then
         write( text,'(a,i4)' )'Main loop at timestep=',na
         call print_memuse_stats(text)
       endif

    enddo

    call mpp_clock_end (id_loop)

!   ------ end of atmospheric time step loop -----

 call atmos_model_end
 call fms_io_exit
 call fms_end

contains

!#######################################################################

   subroutine atmos_model_init

!-----------------------------------------------------------------------
    integer :: unit, ierr, io, logunit
    integer :: ntrace, ntprog, ntdiag, ntfamily
    integer :: date(6)
    type (time_type) :: Run_length
!$    integer :: omp_get_thread_num
    integer :: get_cpu_affinity, base_cpu
    integer :: yr, mo, total_days, total_seconds
!-----------------------------------------------------------------------
!----- initialization timing identifiers ----

 id_init = mpp_clock_id ('MAIN: initialization', grain=CLOCK_COMPONENT)
 id_loop = mpp_clock_id ('MAIN: time loop'     , grain=CLOCK_COMPONENT)
 id_end  = mpp_clock_id ('MAIN: termination'   , grain=CLOCK_COMPONENT)

 logunit = stdlog()

 call mpp_clock_begin (id_init)

!-------------------------------------------
! how many tracers have been registered?
!  (will print number below)
   call register_tracers ( MODEL_ATMOS, ntrace, ntprog, ntdiag, ntfamily )


!----- read namelist -------

#ifdef INTERNAL_FILE_NML
     read (input_nml_file, nml=main_nml, iostat=io)
     ierr = check_nml_error(io, 'main_nml')
#else
   unit = open_namelist_file ( )
   ierr=1; do while (ierr /= 0)
          read  (unit, nml=main_nml, iostat=io, end=10)
          ierr = check_nml_error (io, 'main_nml')
   enddo
10 call mpp_close (unit)
#endif

!----- write namelist to logfile -----

   call write_version_number (version,tag)
   if ( mpp_pe() == mpp_root_pe() ) write (logunit, nml=main_nml)

   if (dt_atmos == 0) then
     call error_mesg ('program atmos_model', 'dt_atmos has not been specified', FATAL)
   endif

   if(lowercase(calendar) == 'no_calendar') then
     calendartype = NO_CALENDAR
   else if(lowercase(calendar) == 'thirty_day_months') then
     calendartype = THIRTY_DAY_MONTHS
   else if(lowercase(calendar) == 'noleap') then
     calendartype = NOLEAP
   else if(lowercase(calendar) == 'julian') then
     calendartype = JULIAN
   else if(lowercase(calendar) == 'gregorian') then
     calendartype = GREGORIAN
   else
     call error_mesg ('program atmos_model', trim(calendar)//' is an invalid value for calendar', FATAL)
   endif
   call set_calendar_type(calendartype)

!----- read restart file -----

   if (file_exist('INPUT/atmos_model.res')) then
       call mpp_open (unit, 'INPUT/atmos_model.res', action=MPP_RDONLY, nohdrs=.true.)
       read  (unit,*) date
       call mpp_close (unit)
   else
    ! use namelist time if restart file does not exist
      if(calendartype == NO_CALENDAR) then
        date(1:2) = 0
        date(3:6) = current_time
      else
        date(1:3) = 1
        date(4:6) = current_time(2:4)
      endif
   endif

!----- write current/initial date actually used to logfile file -----

    if ( mpp_pe() == mpp_root_pe() ) then
      write (logunit,16) date
    endif

 16 format ('  current time used = ',i4,'-',i2,'-',i2,1x,i3,2(':',i2.2)) 

!  print number of tracers to logfile
   if (mpp_pe() == mpp_root_pe()) then
        write (logunit, '(a,i3)') 'Number of tracers =', ntrace
        write (logunit, '(a,i3)') 'Number of prognostic tracers =', ntprog
        write (logunit, '(a,i3)') 'Number of diagnostic tracers =', ntdiag
   endif

!-----------------------------------------------------------------------
!------ initialize diagnostics manager ------

    call diag_manager_init

!----- always override initial/base date with diag_manager value -----

!----- get the base date in the diag_table from the diag_manager ----
!      this base date is typically the starting date for the
!      experiment and is subtracted from the current date

    call get_base_date ( date_init(1), date_init(2), date_init(3), &
                         date_init(4), date_init(5), date_init(6)  )

!----- set initial and current time types ------
!----- set run length and compute ending time -----
#ifdef MARS_GCM
!               Dont allow year, month or minutes in the Mars model
    if ( date_init(1) /= 0 .or. date_init(2) /= 0 .or. date_init(5) /= 0) then
         call error_mesg ('program atmos_model', 'invalid base base - '// &
                          'year, month and minutes must all be zero in the Mars model', FATAL)
    endif
#endif MARS_GCM
    if(calendartype == NO_CALENDAR) then
       Time_init  = set_time(date_init(4)*int(SECONDS_PER_HOUR)+date_init(5)*int(SECONDS_PER_MINUTE)+date_init(6),date_init(3))
       Time       = set_time(date     (4)*int(SECONDS_PER_HOUR)+date     (5)*int(SECONDS_PER_MINUTE)+date     (6),date     (3))
       Run_length = set_time(       hours*int(SECONDS_PER_HOUR)+     minutes*int(SECONDS_PER_MINUTE)+     seconds,days        )
    else
       Time_init  = set_date(date_init(1),date_init(2), date_init(3),date_init(4),date_init(5),date_init(6))
       Time       = set_date(date(1),date(2),date(3),date(4),date(5),date(6))
       Time_tmp = Time
       total_days = 0
       do yr=1,years
         days_tmp = days_in_year(Time_tmp)
         total_days = total_days + days_tmp
         Time_tmp = Time_tmp + set_time(0,days_tmp)
       enddo
       do mo=1,months
         days_tmp = days_in_month(Time_tmp)
         total_days = total_days + days_tmp
         Time_tmp = Time_tmp + set_time(0,days_tmp)
       enddo
       total_days = total_days + days
       total_seconds = hours*3600 + minutes*60 + seconds
       Run_length    = set_time (total_seconds,total_days)
    endif
    Time_end   = Time + Run_length
!-----------------------------------------------------------------------
!----- write time stamps (for start time and end time) ------

      call mpp_open (unit, 'time_stamp.out', form=MPP_ASCII, action=MPP_OVERWR, &
                     access=MPP_SEQUENTIAL, threading=MPP_SINGLE, nohdrs=.true. )

      if ( mpp_pe() == mpp_root_pe() ) write (unit,20) date

!     compute ending time in days,hours,minutes,seconds
      if(calendartype == NO_CALENDAR) then
        call get_time ( Time_end, date(6), date(3) )  ! gets sec,days
        date(4) = date(6)/int(SECONDS_PER_HOUR); date(6) = date(6) - date(4)*int(SECONDS_PER_HOUR)
#ifdef MARS_GCM
        date(5) = 0                                ; date(6) = date(6) - date(5)*int(SECONDS_PER_MINUTE)
#else
        date(5) = date(6)/int(SECONDS_PER_MINUTE)  ; date(6) = date(6) - date(5)*int(SECONDS_PER_MINUTE)
#endif
      else
        call get_date(Time_end, date(1), date(2), date(3), date(4), date(5), date(6))
      endif
      if ( mpp_pe() == mpp_root_pe() ) write (unit,20) date

      call mpp_close (unit)

  20  format (6i7,2x,'day')   ! can handle day <= 999999

!-----------------------------------------------------------------------
!--- compute the time steps ---
!    determine number of iterations through the time integration loop 
!    must be evenly divisible

      Time_step_atmos = set_time (dt_atmos,0)
      num_atmos_calls = Run_length / Time_step_atmos

!-----------------------------------------------------------------------
!----- initial (base) time must not be greater than current time -----

   if ( Time_init > Time ) call error_mesg ('program atmos_model',  &
                   'initial time is greater than current time', FATAL)

!----- make sure run length is a multiple of atmos time step ------

   if ( num_atmos_calls * Time_step_atmos /= Run_length )  &
        call error_mesg ('program atmos_model',  &
           'run length must be multiple of atmosphere time step', FATAL)
   
!-----------------------------------------------------------------------
!------ initialize atmospheric model ------

!$      call omp_set_num_threads(atmos_nthreads)
      if (mpp_pe() .eq. mpp_root_pe()) then
        unit=stdout()
        write(unit,*) ' starting ',atmos_nthreads,' OpenMP threads per MPI-task'
        call flush(unit)
      endif
      base_cpu = get_cpu_affinity()
!$OMP PARALLEL
!$      call set_cpu_affinity(base_cpu + omp_get_thread_num())
#ifdef DEBUG
!$      write(6,*) 'PE: ',mpp_pe(),'  thread_num', omp_get_thread_num(),'  affinity:',get_cpu_affinity()
!$      call flush(6) 
#endif
!$OMP END PARALLEL

      call atmosphere_init (Time_init, Time, Time_step_atmos)
      call atmosphere_domain(atmos_domain)

!-----------------------------------------------------------------------
!   open and close dummy file in restart dir to check if dir exists
      call mpp_set_current_pelist()
      call mpp_open  (unit, 'RESTART/file' )
      call mpp_close (unit, action=MPP_DELETE)

!  ---- terminate timing ----
   call mpp_clock_end (id_init)

!-----------------------------------------------------------------------

   call print_memuse_stats('atmos_model_init')
   end subroutine atmos_model_init

!#######################################################################

   subroutine atmos_model_end

   integer :: unit, date(6)
!-----------------------------------------------------------------------
   call mpp_clock_begin (id_end)

      call atmosphere_end

!----- compute current time in days,hours,minutes,seconds -----

      if(calendartype == NO_CALENDAR) then
        date(1:2) = 0
        call get_time ( Time, date(6), date(3) )
        date(4) = date(6)/int(SECONDS_PER_HOUR); date(6) = date(6) - date(4)*int(SECONDS_PER_HOUR)
#ifdef MARS_GCM
        date(5) = 0                              ; date(6) = date(6) - date(5)*int(SECONDS_PER_MINUTE)
#else
        date(5) = date(6)/int(SECONDS_PER_MINUTE); date(6) = date(6) - date(5)*int(SECONDS_PER_MINUTE)
#endif MARS_GCM
      else
        call get_date(Time, date(1), date(2), date(3), date(4), date(5), date(6))
      endif

!----- check time versus expected ending time ----

      if (Time /= Time_end) call error_mesg ('program atmos_model',  &
              'final time does not match expected ending time', WARNING)

!----- write restart file ------

      if ( mpp_pe() == mpp_root_pe() ) then
           call mpp_open (unit, 'RESTART/atmos_model.res', form=MPP_ASCII, action=MPP_OVERWR, &
                          access=MPP_SEQUENTIAL, threading=MPP_SINGLE, nohdrs=.true. )
           write (unit,'(6i6,8x,a)') date, &
                 'Current model time: year, month, day, hour, minute, second'
           call mpp_close (unit)
      endif

!----- final output of diagnostic fields ----
      call set_domain(atmos_domain)  ! This assumes all output fields are on the atmos domain

      call diag_manager_end (Time)

      call nullify_domain()

      call mpp_clock_end (id_end)
!-----------------------------------------------------------------------

   end subroutine atmos_model_end

!#######################################################################
! routines to set/get date when no calendar is set (i.e., yr=0 and mo=0)
!#######################################################################

end program atmos_model
